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When DNA molecules are heated they denature. This occurs locally so that loops of molten 
single DNA strands form, connected by intact double-stranded DNA pieces. The properties of this 
‘melting’ transition have been intensively investigated. Recently there has been a surge of interest in 
this question, caused by experiments determining the properties of partially bound DNA confined to 
nanochannels. But how does such confinement affect the melting transition? To answer this question 
we introduce, and solve a model predicting how confinement affects the melting transition for a 
simple model system by first disregarding the effect of self-avoidance. We find that the transition 
is smoother for narrower channels. By means of Monte-Carlo simulations we then show that a 
model incorporating self-avoidance shows qualitatively the same behaviour and that the effect of 


confinement is stronger than in the ideal case. 

PACS numbers: 87.15.A-,87.15.Fh,36.20.Ey,87.14.gk 

Thermal denaturation of DNA molecules has been 
studied intensively in the past decades [1-6]. This pro¬ 
cess, also referred to as ‘melting’, describes how the two 
strands of the double-stranded helix separate when the 
hydrogen bonds holding together the base pairs disinte¬ 
grate. When this happens at not too high temperatures, 
loops or DNA ‘bubbles’ of different sizes form. The fun¬ 
damental question is to determine which factors, apart 
from temperature, influence the shape of the loop-size 
distribution, and how this distribution varies as the tem¬ 
perature is changed. A first step is to analyse a sim¬ 
ple summary statistic - the average fraction of unbroken 
bonds. This binding probability is the order parameter 
of the melting transition. 

There are two important factors determining the tem¬ 
perature dependence of this order parameter. First, G-C 
bonds are stronger than A-T bonds. Since A-T-rich re¬ 
gions therefore melt more easily [4], DNA sequence influ¬ 
ences the melting transition. Second, entropy plays an 
important role. Loops and bound DNA strands differ in 
how the number of spatial configurations depends on the 
number of base pairs, and this difference determines the 
properties of the melting transition. 

Recently there has been great interest in DNA melt¬ 
ing, stimulated by the possibility to optically study 
sequence-specific local melting of DNA molecules con¬ 
fined to nanochannels, see Ref. [7] for a review. The 
aim of these studies is to experimentally distinguish G- 
C-rich regions from A-T-rich regions, obtaining sequence- 
dependent ‘bar codes’ representing given DNA sequences. 
Potential applications of this method include fast species 
and strain identification [7]. In the experiments described 
in Refs. [7, 8] DNA is confined to a channel that is signif¬ 
icantly smaller than the radius of gyration of the uncon¬ 
fined DNA molecule. In this case confinement substan¬ 
tially affects the conformational fluctuations of the DNA 
molecule. This indicates that confinement must have a 


strong effect upon the melting transition. For very short 
DNA molecules (15 base-pairs), the effect of confining the 
molecules to channels with a size of a few nanometers has 
recently been studied by molecular dynamics simulations 
[9]. Yet for longer molecules which exhibit loops of many 
different sizes, it is not known how confinement influences 
the melting transition. 

This motivated us to analyse an idealised model system 
describing the melting of a long double-stranded chain 
confined to a channel of width D. In its simplest form 
the model disregards self-avoidance. In this case we can 
quantitatively explain the effect of confinement upon the 
melting transition. We find that the temperature de¬ 
pendence of the order parameter is smoother for more 
strongly confined chains. Our theory explains this ef¬ 
fect by considering how confinement affects the config¬ 
urational entropy of single-stranded loops. By means 
of Monte-Carlo simulations we show that self-avoiding 
chains show a similar trend: confinement renders the 
melting transition smoother. The effect is much stronger 
than in the ideal case (disregarding self-avoidance). 

The statistics of melting can be described by treating 
the DNA molecule as a sequence of molten and intact 
sections. This method is often referred to as the Poland- 
Scheraga model [Ij. Within this model, the statistics of 
DNA melting depend on the entropy of a single-stranded 
loop. This entropy is commonly assumed to be deter¬ 
mined by the return probability of a random walk repre¬ 
senting a single DNA strand closing upon itself [1, 10-12]. 
Our analysis shows that this is, in general, incorrect. We 
demonstrate that the relevant quantity is instead the first 
return probability [13]. For unconfined random walks in 
three dimensions this makes little difference since these 
probabilities have the same power-law dependence upon 
loop size [13]. Yet it is crucial for determining the en¬ 
tropy of large, confined loops. The entropy of large loops 
in turn determines the order of the phase transition. By 
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FIG. 1. (a) Conformation of a partially molten DNA molecule confined to a channel of width D — 14^k. Shown is a snapshot 

obtained from computer simulations of the model described in the text. The parameters of the simulation are N = 400, r = Iyl, 
a = 0.63^k, Fib = 2.4. The beads of diameter a are not drawn. Colour code: Molten parts (loops) are red, Intact parts are blue, 
(b) Schematic illustration of the three-dimensional model of a confined double- stranded chain that is used in the Monte-Carlo 
simulations, see text. Note that the two strands are clamped at the ends. 


analysing this entropy, we show that the order of the 
phase transition is unchanged by confinement in a chan¬ 
nel. This result is consistent with a study by Causo et al. 
[14], who show that for a lattice model of DNA melting, 
the order of the phase transition is the same in one and 
three dimensions. Yet as we discuss above, the melting 
transition is strongly influenced by confinement. We in¬ 
fer that the order of the phase transition describes only a 
minute range of temperatures, whereas the overall shape 
of the temperature dependence of the order parameter 
depends strongly on the value of first return probability 
for smaller loops. 

Model. We treat a highly idealised model of double- 
stranded DNA, illustrated in Fig. 1. Each strand consists 
of a freely jointed chain of N spherical beads of diameter 
a and excluded volume v = 47ra^/3, connected by ideal 
rods of length £k- The monomers of the first chain are 
labeled by n = 1, 2,..., A^. The monomers of the second 
chain are labeled by n' (same range). The two chains 
are attached to each other at both ends, as shown in 
Fig. 1(b). For a given monomer n and its partner n' = n 
(for 1 < n = n' < N) we test whether the distance 
between the two monomers is less than r. If so the two 
monomers are considered bound, and a negative binding 
energy — Ab (in units such that k^T = 1) is assigned 
to this pair. Note that binding is only allowed between 
monomers corresponding to different chains, and only for 
n = n'. 

Effect of confinement on melting. The statistics of 
melting is determined by the difference in free energy 
between a bound section of m base pairs and a molten 
section of the same size, bracketed by closed base pairs 
[11]. The bound section has a higher energy, but also a 
higher entropy, and the probability of melting is deter¬ 
mined by the interplay of the two. Within our model the 
increase in energy upon melting is given by the number 
of broken bonds between the bases, proportional to m. 
The increase in entropy, on the other hand, has a more 


complicated dependence upon m. To figure out this de¬ 
pendence it is helpful to consider the two chains of the 
DNA as independent chains, joined at one end. For the 
chains to remain bound they must move in lockstep. The 
molten state is less restrictive, allowing any configuration 
of the two strands that fulfill two conditions. First, the 
two strands must together from a closed loop, thus they 
must end at the same position. Second, the loop must 
not close prematurely, as that would not correspond to 
a single molten region of size m, but rather two or more 
smaller molten regions. The second condition has not 
been considered in previous studies of melting [1, 10-12]. 
From this condition it follows that the entropy of a loop 
is not given by the return probability of a random walk, 
but by its first return probability. 

The first step in deriving a theory for melting is to 
determine the melting entropy AS{m)/kB by counting 
the number of configurations Db and flu satisfying the 
conditions of the bound and molten state, respectively. 
For the ideal case of our model (a = 0, no self-avoidance) 
we proceed as follows. For the bound state, we estimate 
Db « Here c is the number of ways in which 

a single Kuhn length segment can be oriented. For a 
spatially continuous model such as ours, it is in principle 
infinite. Imagine that space is finely discretised. Then 
c is finite and we can show that c drops out in the final 
result. Furtherp = P(jrfc —r^j < r and < r) 

is the probability that monomer k remains in lockstep, 
given that monomer F — 1 is in lockstep, and Kim) is the 
probability that a chain of length m does not leave the 
channel. This probability can be determined from the 
solution of a diffusion equation [15, 16]: 

A'(to) « exp [—7r^F|-m/(3Il^)] . (1) 

For the molten state, we find Du K(m)^ f (m). The 

function /(to) is the probability that the two chains form 
a single closed loop, or alternatively the probability that 
the random walk performed by the separation r^ —r). first 
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FIG. 2. A sketch of the first return probability f{m). Dashed 
line: f{m) for an unconfined DNA molecule. Solid line: f{m) 
for ideal DNA confined to a square channel of size D ^ Ik- 

returns to within a radius r of the origin after m steps. 
The function f{m) is sketched in Fig. 2. For small values 
of m the probability f{m) is simply the first return prob¬ 
ability for an unconfined three-dimensional random walk 
- a short loop does not feel the presence of the walls. As¬ 
suming that m ^ 1, the first return probability scales as 
/(to) oc to“^/^ in this region [13]. In the second and third 
regions the problem is essentially one-dimensional. Con¬ 
sider the component of the random walk in the channel 
direction. This component performs a one-dimensional 
random walk. Each time it returns to the origin, there is 
a probability k oc /{D'^iK) that the three -dimensional 
random walk returns to within a radius r of the origin 
[16, 17]. The probability of first return can therefore be 
mapped to the solved problem of finding the probability 
of absorption for a one-dimensional random walk with a 
partially absorbing sink at the origin [18]. One finds that 
this probability scales as /(to) oc for small to, and 

as /(to) oc for large to. In summary, the func¬ 

tion /(to) has the form sketched in Fig. 2. In the limit 
r —> 0 , one need not distinguish between the probability 
of return and first return. We discuss this limit in the 
supplementary material. 

Taking the above results together we have 

O 77? 

Ag/fcB = log ^ = a-byTO- -blog/( to). (2) 

Here a. and 7 are undetermined constants which do not 
depend on whether the DNA is confined or not, a de¬ 
pends on the undetermined prefactors in our expressions 
for Hb and Hu, and 7 = log(l/p) depends on the probabil¬ 
ity for a bound section to stay in lockstep. The terms that 
are proportional to to have the same effect as shifting the 
binding energy according to e = Ab — 7 -I- tt^ 4to/(3D2). 

Given Eq. (2) we can compute the order parameter, 
the probability pb that a base pair is bound, at different 
binding energies. Let us first consider the qualitative be¬ 
haviour of the melting curve in the thermodynamic limit 
[N —>■ 00 ). It is well known that the melting probability 
can be computed from a grand-canonical description (see 


for example Refs. [1, 10, 11]), 


= [l + e^+'^*z;(p*)]-i, 

(3a) 

00 

= e“^e”^^*/(™), 

(3b) 

771=1 

= e^+'^'[l + Z,(p*)]. 

(3c) 


Here is the grand-canonical partition function of 

a molten section with chemical potential p, and Z'^ is 
its derivative with respect to p. Eq. (3c) specifies the 
value of the chemical potential necessary to reach the 
thermodynamic limit. If Eq. (3c) admits no solution for 
real p*, we have pb = 0 . 

The effect of confinement upon the order parameter is 
determined by /(to). The D-dependent correction to e 
is small because we assume that Ik ^ D. At high values 
of e, Eq. (3c) yields /i* « —e ^ 0. For such large nega¬ 
tive values of p*, the sum in Eq. (3b) converges rapidly. 
The physical interpretation is that the formation of large 
molten loops is very unlikely when is large. The or¬ 
der parameter is therefore determined by the behaviour 
of /(to) for small values of to. Yet here /(to) is not in¬ 
fluenced by confinement as Fig. 2 shows. In this region, 
then, confinement influences the melting curve only very 
slightly. In the opposite limit of negative values of e, 
the fact that /(to) is significantly larger for the confined 
chain at large to implies that Eq. (3c) admits solutions 
for a larger range of binding energies, and thus that pb 
increases at lower energies for the confined chain, com¬ 
pared to the unconhned one. At intermediate binding 
energies, the binding probability depends sensitively on 
the exact shape of /(to), and is therefore hard to predict. 
But we can conclude that the melting curve is sharper for 
the unconhned chain, and smoother for the confined one. 

Note that although the shape of the melting curve is 
smoother for the confined chain, the order of the phase 
transition must remain unchanged, as the asymptotic 
scaling (/(to) oc to“^/^ as to — >■ 00 ) is the same for 
the confined and the unconfined chain. This shows that 
the order of the phase transition only determines the be¬ 
haviour of the melting curve for a very small range of 
binding energies close to the point where pb becomes non¬ 
zero. Yet determining the location of this point and the 
full shape of the melting curve requires one to consider 
the behaviour of /(to) for all to. 

Simulation results. We performed computer simula¬ 
tions of our model using the Metropolis algorithm with 
crankshaft trial updates [19]. A snapshot from such sim¬ 
ulations is shown in Fig. 1(a). We begin by comparing 
the results for the ideal case (a = 0 ) to the theoretical 
predictions above. Fig. 3(a) shows the probability Ph{n) 
that a given base pair is closed, for an ideal chain with 
binding energy — 0.78, at different levels of confine¬ 
ment. For the upper two curves, this binding probability 
is independent of n except close to the ends of the chain. 
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FIG. 3. Effect of confinement upon the melting transition, 
ideal case, (a) Binding probability as a function of monomer 
number, for an ideal chain with binding energy Eb = 0.78, at 
different levels of confinement: D/^k = 8 (blue), D/£k = 14 
(red), unconfined (green). N = 400, a = 0, r = £k (b) Order 
parameter ph{N/2) as a function of the binding energy, for 
different channel sizes as in panel (a). Solid lines: Fits to a 
solution of Eqs. (3), see main text for details. 


This indicates that finite size effects are not important 
for these binding energies. The lower curve shows a case 
where the binding probability does not have a clearly 
plateau in the middle. In this case, then, finite size ef¬ 
fects must result in a measured binding probability that 
is larger than it would be for an infinite chain. This effect 
is only present for small binding energies, and is largest 
for the unconfined chain. The simulations thus underes¬ 
timate the sharpness of the unconfined melting curve. 

Fig. 3(b) shows how confinement influences the melt¬ 
ing transition. We plot the probability that the middle 
base pair is closed, as a function of the binding energy 
E\y. We find that confinement does not influence the 
melting probability at large i?bi but also that increas¬ 
ing confinement leads to a smoother transition overall. 
These results are in perfect agreement with the theoret¬ 
ical expectations discussed above. 

We also show that the melting curve of Fig. 3(b) can 
be well described by a numerical solution to Eq. (3). 
For this numerical solution, we assume that f{m) = 
(m -I- consistent with the shape 

shown in Fig. 2. The empirical parameter A is required to 
compensate for the fact that the diffusion approximation 
fails for very short loops. For m > the func¬ 
tion f{m) crosses over to the scaling fim) oc A 

is another fitting parameter, which reflects the fact that 
the exact location of the cross-over between and 

is not known. Since the polymers that we have 
simulated are too short to exhibit the third region, we do 



n 



FIG. 4. Effect of confinement upon the melting transition, 
self-avoiding case, (a) Binding probability as a function of 
monomer number, for a self-avoiding chain with binding en¬ 
ergy Eb = 2.33, at different levels of confinement: DjlK = 8 
(blue), DI£k = 14 (red), unconfined (green). N = 400, 
a = 0.63^k, r = (b) Order parameter pb(lV/2) as a func¬ 

tion of the binding energy, for different channel sizes as in 
panel (a). 

not include it in our numerical model. Using this expres¬ 
sion for /(m), Eqs. (3) can be evaluated in Mathemat- 
ica, yielding closed-form expressions forpb(M*) and e(p*). 
Varying fi* from —oo to 0 then allows one to trace out 
the curve Pb(e)- The parameters a. A, A, and 7 were de¬ 
termined using the Mathematica routine FindFit. Eirst, 
the values a = —0.24, 7 = 1.58 and A = 0.97 were found 
by fitting the unconfined curve. Then, keeping these val¬ 
ues fixed, the value A = 2.99 was found by fitting the 
curve corresponding io D = 8 £k- Although there are 
four fitting parameters, the fact that the same parameter 
values yield good fits at different confinement strengths 
confirms that our analysis accurately describes the effect 
of confinement on the melting transition. That the over¬ 
lap is not exact is to be expected, as the function /(m) is 
only approximately correct. Note also that whereas the 
theoretical curve for the unconfined chain shows a kink 
at low bending energies, the simulated curve is smooth. 
This discrepancy is caused by the finite size of the simu¬ 
lated system (see discussion above). 

Effect of self avoidance. Eor the case of self-avoiding 
chains we do not have theoretical predictions. However, 
we have performed simulations of the model described 
above with a = 0.63£k- The results are shown in Fig. 4. 
Fig. 4(b) shows that confinement has the same qualitative 
effect here as for the ideal chain, in that it makes the tran¬ 
sition smoother. Yet the effect is even more noticeable 
in this case. Note that the for the case of the unconfined 
chain (green curve), the transition is sharper than for 
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an ideal chain [see Fig. 3(b)]. This is in agreement with 
established theory for unconfined chains [10, 11]. The 
fact that the lower curve in Fig. 3(a) does not exhibit 
a clear plateau indicates that for the unconfined chain, 
finite size effects influence the melting probability at the 
transition. Also, compared to the ideal case, much larger 
binding energies are required for binding to occur. This 
is a consequence of the fact that it is much harder for a 
bound region of the self-avoiding chain to remain in lock- 
step, as a large part of the binding volume is inaccessible 
because of self-avoidance. 

Discussion. We have studied how the melting transi¬ 
tion of DNA is modified by confinement. To this end, 
we have analysed a simple model by analytical calcula¬ 
tions and Monte Carlo simulations. We show that con¬ 
finement has a profound effect on the entropy of large 
molten sections, and that this results in a smoother melt¬ 
ing transition than in the unconfined case. We show that 
the melting transition can be successfully analysed within 
the Poland-Scheraga model by calculating the entropy of 
a confined loop, and that this entropy is determined by 
the first return probability of a random walk. Our theory 
is in quantitative agreement with simulation results for 
an ideal polymer model. Simulations of a self-avoiding 
polymer show a similar, but even stronger effect of con¬ 
finement. 

The model that we study is clearly not a very realis¬ 
tic description of DNA. Yet the effect of confinement on 
the entropy of molten sections is insensitive to the details 
of the model. We therefore expect that channel confine¬ 
ment has a similar qualitative effect upon the melting 
transition, i.e. that the transition is smoother for con¬ 
fined DNA. The main difference is that for real DNA, the 
stiffness of dsDNA is much higher than that of ssDNA. 
Thus K{m) [Eq. (1)] is different for a bound section than 
for a molten one. This has the effect of shifting the melt¬ 
ing transition to lower temperatures, i.e. in the opposite 
direction compared to our model. We explore the conse¬ 
quences of our theory for real DNA in Ref. [20]. 
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In this supplementary material, we derive an explicit 
formula for the difference in entropy AS'(m) between a 
bound and a molten region, in the limit r —)■ 0. We show 
that it agrees with the scaling in Eq. (2) of the main text, 
where /(m) obeys the scaling of the first two regions of 
Fig. 2 (the third region, where for large m f{m) oc rnT^I'^, 
does not appear when r —>■ 0.) We further note that while 
the exact expression of the formula depends on how the 
region is assumed to be connected to the rest of the DNA, 
the scaling of Eq. (2) holds for all cases. 

Denote by Gjv(r; ro)d^r the probability that a chain 
started at rg survives the addition of JV monomers and 
ends up in the volume d^r around r. If the channel size is 
much larger than the Kuhn length of the chain, G'Ar(r; rg) 
obeys a diffusion equation [ 1 , 2 ] 


dGjv(r; rg) 
dN 


p2 


VrGAr(r;ro), 


with absorbing boundary conditions: GAr(r;rg) = 0 at 
the walls. If the channel is square, with side length D, 
the solution can be found by separation of variables [I, 2] 


G 7 v(r;ro) = G^{x; Xo)G^{y;yo)G^l{z; zq), where 
G%{z;zo) = ^j3/{2TTNel)exp ) . 


kTTXo 


knx 


Gn{x; xo) = 2^ sm —^ sm — exp ^ - 


k=l 


D 


D 


6 D 2 


Since the 2 ;-direction is unconstrained, G^l^{z;zo) is sim¬ 
ply the probability distribution of a one-dimensional ran¬ 
dom walk. 


From these expressions we obtain the fraction of all 
chains which “survive” for N links by integrating over 
all start and end positions: 

Gn = JJ ^^j^d^rGN{r;ro) (SI) 

2 f eiTrH‘2k + l)^N\Y 

(2fc + I)7rJ 6D^ ) ■ 

This expression is approximately proportional to the 
number of configurations of a bound segment under con¬ 
finement, Db(m) « c^'^p^Gra- See the main text for a 
discussion of the parameters c, p. 




Let us now turn to the other quantity of interest: 
Du(m), the number of possible configurations for a con¬ 
fined loop. The probability that a chain started at r 
survives the addition of N links and returns to within 
d^r of the starting point is given by 


GAr(r; r)d^r = G^ix; x)Gji{y; y)G%{z; z)d^v 


= y 3 /( 27 rIV 4 ) 


E 2 . 2 klTX I 

sm ^exp - 


.fc=i 


/ 


E 


2 . 2 kiry 

— sm -exp 1 — 

D D ^ ' 


6 D 2 


d-^r. 


Averaging the return probability over all starting posi¬ 
tions r yields 


G% = ^5/i27rml)D-^ 


oo 


Eexp 

.fc=i 




2 

• (S2) 


This quantity is proportional to the number of configura¬ 
tions for a ring-polymer confined to a square nanochannel 
of diameter D, Du(m) oc 

The entropy difference between a bound and a molten 
region is given by AS{m)/kB = log(Du/Db). We now 
show that this expression for AS agrees with Eq. (2) of 
the main text, with /(m) scaling as the first two regions 
of Fig. 2 (recall that the third region does not occur when 
r —> 0, which is the limit we treat here). Extracting 
the dominant exponential terms from Eqs. (SI) and (S2) 
yields 


^ 2 /i 2 , 


AS/kB = a + ym — 
-blog 


3D2 
Du( 2 m) exp 


log 


Db(m)exp 


2£^7r^fc^m\ 
3D2 ) 

^7r^A:^m\ 
3D2 ) ■ 


(S3) 


As in Eq. (2) of the main text, a and 7 are undetermined 
constants. The last term in the expression is approxi¬ 
mately constant, the argument of the logarithm stays in 
the interval [1-1.53] for all m. It remains to consider 
the expression F{m) = Du(2m) exp[ 2 £|^ 7 r^fc^m/( 3 D^)]. 
For large m » only the prefactor in Eq. (S2) re¬ 
mains, yielding F{m) oc In the opposite limit 

m Poisson resummation shows that F{m) oc 
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m^-‘3/2^ In summary, the function F{m) exhibits the same 
scaling as /(to) in the limit r —)• 0, and thus Eq. (S3) re¬ 
produces exactly Eq. (2) of the main text. 

In the derivation of Eq. (S3), we averaged Gn over 
different positions, with all positions in the channel 
weighted equally. While this procedure gives the correct 
statistics for an isolated polymer, it is not quite correct 
for a region in the interior of a larger polymer. Instead, 
the averaging must give more weight to positions which 
are far from the channel walls [2]. How the weighting is to 
be done depends on the state of the surrounding polymer. 
If the region in question is in the interior of an otherwise 
linear polymer the end position has a weight proportional 
to sin(7ra:/D) sin(7ry/£)) [2]. If, on the other hand, the 
region is bracketed by large molten regions (loops), the 
weight is instead proportional to sin^(7rx/Z)) sin^(7rj//D). 
It thus appears that deriving an exact formula for AS'(to) 
is not possible. However all three averaging schemes dis¬ 
cussed here show the same general behaviour, given by 
Eq. (2) and Fig. 2 of the main text. 


Note that in the limit r —)■ 0 the probability of return 
and the probability of first return are equal. Therefore, 
the third region in Fig. 2 of the main text cannot be ob¬ 
tained from the calculations described here. However, for 
small but non-zero values of r it is easy to see that when¬ 
ever the z-coordinates of two complementary monomers 
are close then there is a small, finite probability for bind¬ 
ing to occur. As stated in the main text, this allows 
one to map the problem of calculating the probability of 
first return to a solved one-dimensional problem [3], from 
which the existence of the third region in Fig. 2 of the 
main text follows. 
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